library(foreign)
library(boot)
CPS <- read.dta("rdcps.dta")
CPS <- subset(CPS, year != 1998 & year != 2000)

################################################################################
## Prepare data
################################################################################
# Dichotomize instrument
CPS$ca_dich <- ifelse(CPS$ca >= 10, 1, 0)

# Control variables
CPS$black <- ifelse(CPS$race==2, 1, 0)
CPS$female <- 1 - CPS$male
CPS$fe_year <- as.factor(CPS$year)
CPS$fe_sob <- as.factor(CPS$birthpl)

################################################################################
## Model formulae
################################################################################
eq1 <- hs ~ ca_dich + black + female + fe_sob + age + age2 + age3 + age4
eq2 <- vote ~ hs + black + female + fe_sob + age + age2 + age3 + age4
eq3 <- vote ~ hshat + black + female + fe_sob + age + age2 + age3 + age4
inst <- ~ ca_dich + black + female + fe_sob + age + age2 + age3 + age4

################################################################################
## Bootstrap function
################################################################################
paboot.2sls <- function(data, index){
    # Subset data for bootstrap
    D <- data[index,]
    
    # Fit 2sls
    fit1 <- lm(eq1, D)
    D$hshat <- fit1$fitted
    fit2 <- lm(eq3, D)

    # Predicted probabilities
    Z1 <- Z0 <- fit1$model
    Z1[,"ca_dich"] <- 1
    Z0[,"ca_dich"] <- 0

    PX1 <- pmax(0, pmin(1, predict(fit1, newdata=Z1)))
    PX0 <- pmax(0, pmin(1, predict(fit1, newdata=Z0)))

    X1 <- X0 <- data.frame(fit2$model)
    X1[,"hshat"] <- 1
    X0[,"hshat"] <- 0

    PY1 <- pmax(0, pmin(1, predict(fit2, newdata=X1)))
    PY0 <- pmax(0, pmin(1, predict(fit2, newdata=X0)))

    P11.1 <- PY1 * PX1
    P10.1 <- PY0 * (1 - PX1)
    P01.1 <- (1 - PY1) * PX1
    P00.1 <- (1 - PY0) * (1 - PX1)
    P11.0 <- PY1 * PX0
    P10.0 <- PY0 * (1 - PX0)
    P01.0 <- (1 - PY1) * PX0
    P00.0 <- (1 - PY0) * (1 - PX0)

    Q <- mean(fit1$model$ca_dich)

    P11 <- P11.1 * Q + P11.0 * (1-Q)
    P10 <- P10.1 * Q + P10.0 * (1-Q)
    P01 <- P01.1 * Q + P01.0 * (1-Q)
    P00 <- P00.1 * Q + P00.0 * (1-Q)

    # Compute PCA & CPCA estimates
    ITT <- P11.1 + P10.1 - P11.0 - P10.0

    pca.MonoXY <- mean(pmax(0, pmin(1, ITT/(P11.1 - P11.0), na.rm=T)))
    pca.MonoX.lo <- mean(pmax(0, pmin(1, ITT/(P11.1 - P11.0), na.rm=T)))
    pca.MonoX.up <- mean(pmin(1, pmax(0, (P00.0 - P00.1)/(P11.1 - P11.0), na.rm=T)))
    pa.MonoXY.lo <- mean(pmax(0, pmin(1, ITT * Q / P11, na.rm=T)))
    pa.MonoXY.up <- mean(pmin(1, pmax(0, (ITT * Q + P11.0) / P11, na.rm=T)))
    pa.MonoX.lo <- mean(pmax(0, pmin(1, ITT * Q / P11, na.rm=T)))
    pa.MonoX.up <- 1
    late <- coef(fit2)["hshat"]
    c(pca.MonoXY, pca.MonoX.lo, pca.MonoX.up, 
        pa.MonoXY.lo, pa.MonoXY.up, pa.MonoX.lo, pa.MonoX.up, late, Q)
}

################################################################################
## Estimate
################################################################################
res.pooled.full <- boot(CPS, paboot.2sls, 2000)

# Flip 0 and 1 to change QoI
CPS$ca_dich <- 1 - CPS$ca_dich
CPS$hs <- 1 - CPS$hs
CPS$vote <- 1 - CPS$vote

res.pooled.flip <- boot(CPS, paboot.2sls, 2000)

################################################################################
## Plotting
################################################################################
source("IMCI.R")

pdf("fig3-left.pdf", width=4.5, height=4.5)
par(mfrow=c(1,1), oma=c(1,1,1,1), mar=c(3.1,4,2.1,1.1))

point.pca <- res.pooled.full$t0[1]
bounds.pca <- res.pooled.full$t0[2:3]
bounds1.pa <- res.pooled.full$t0[4:5]
bounds2.pa <- res.pooled.full$t0[6:7]

sd.point.pca <- sd(na.omit(res.pooled.full$t[,1]))
sd.bounds.pca <- apply(na.omit(res.pooled.full$t[,2:3]), 2, sd)
sd.bounds1.pa <- apply(na.omit(res.pooled.full$t[,4:5]), 2, sd)
sd.bounds2.pa <- apply(na.omit(res.pooled.full$t[,6:7]), 2, sd)

CI.point.pca <- c(min(1, max(0, point.pca - qnorm(.975)*sd.point.pca)), 
    min(1, max(0, point.pca + qnorm(.975)*sd.point.pca)))
CI.bounds.pca <- pmin(1, pmax(0, IMCI(l=bounds.pca[1], u=bounds.pca[2], 
    var.l=sd.bounds.pca[1]^2, var.u=sd.bounds.pca[2]^2)$ci))
CI.bounds1.pa <- pmin(1, pmax(0, IMCI(l=bounds1.pa[1], u=bounds1.pa[2], 
    var.l=sd.bounds1.pa[1]^2, var.u=sd.bounds1.pa[2]^2)$ci))
CI.bounds2.pa <- pmin(1, pmax(0, IMCI(l=bounds2.pa[1], u=bounds2.pa[2], 
    var.l=sd.bounds2.pa[1]^2, var.u=sd.bounds2.pa[2]^2)$ci))

plot(c(-.1,-.1), bounds1.pa, ylim=c(0,1), xlim=c(-.5,1.5), lwd=5, type="l",
    ylab = "Probability of Causal Attribution", xlab = "",
    main = "Educated Voters Who \n Would Not Have Voted If Uneducated", xaxt="n")
lines(c(-.1,-.1), CI.bounds1.pa)
points(c(-.1,-.1), CI.bounds1.pa, pch="_", cex=1)
lines(c(.1,.1), bounds2.pa, lwd=5)
lines(c(.1,.1), CI.bounds2.pa)
points(c(.1,.1), CI.bounds2.pa, pch="_", cex=1)
mtext("Population", 1, line=1, at=0)

points(.9, point.pca, pch=16)
lines(c(.9,.9), CI.point.pca)
points(c(.9,.9), CI.point.pca, pch="_", cex=1)
lines(c(1.1,1.1), bounds.pca, lwd=5)
lines(c(1.1,1.1), CI.bounds.pca)
points(c(1.1,1.1), CI.bounds.pca, pch="_", cex=1)
mtext("Complier", 1, line=1, at=1)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

dev.off()
rm(res.pooled.full)

################################################################################

pdf("fig3-middle.pdf", width=4.5, height=4.5)
par(mfrow=c(1,1), oma=c(1,1,1,1), mar=c(3.1,4,2.1,1.1))

point.pca <- res.pooled.flip$t0[1]
bounds.pca <- res.pooled.flip$t0[2:3]
bounds1.pa <- res.pooled.flip$t0[4:5]
bounds2.pa <- res.pooled.flip$t0[6:7]

sd.point.pca <- sd(na.omit(res.pooled.flip$t[,1]))
sd.bounds.pca <- apply(na.omit(res.pooled.flip$t[,2:3]), 2, sd)
sd.bounds1.pa <- apply(na.omit(res.pooled.flip$t[,4:5]), 2, sd)
sd.bounds2.pa <- apply(na.omit(res.pooled.flip$t[,6:7]), 2, sd)

CI.point.pca <- c(min(1, max(0, point.pca - qnorm(.975)*sd.point.pca)), 
    min(1, max(0, point.pca + qnorm(.975)*sd.point.pca)))
CI.bounds.pca <- pmin(1, pmax(0, IMCI(l=bounds.pca[1], u=bounds.pca[2], 
    var.l=sd.bounds.pca[1]^2, var.u=sd.bounds.pca[2]^2)$ci))
CI.bounds1.pa <- pmin(1, pmax(0, IMCI(l=bounds1.pa[1], u=bounds1.pa[2], 
    var.l=sd.bounds1.pa[1]^2, var.u=sd.bounds1.pa[2]^2)$ci))
CI.bounds2.pa <- pmin(1, pmax(0, IMCI(l=bounds2.pa[1], u=bounds2.pa[2], 
    var.l=sd.bounds2.pa[1]^2, var.u=sd.bounds2.pa[2]^2)$ci))

plot(c(-.1,-.1), bounds1.pa, ylim=c(0,1), xlim=c(-.5,1.5), lwd=5, type="l",
    ylab = "Probability of Causal Attribution", xlab = "",
    main = "Uneducated Nonvoters Who \n Would Have Voted If Educated", xaxt="n")
lines(c(-.1,-.1), CI.bounds1.pa)
points(c(-.1,-.1), CI.bounds1.pa, pch="_", cex=1)
lines(c(.1,.1), bounds2.pa, lwd=5)
lines(c(.1,.1), CI.bounds2.pa)
points(c(.1,.1), CI.bounds2.pa, pch="_", cex=1)
mtext("Population", 1, line=1, at=0)

points(.9, point.pca, pch=16)
lines(c(.9,.9), CI.point.pca)
points(c(.9,.9), CI.point.pca, pch="_", cex=1)
lines(c(1.1,1.1), bounds.pca, lwd=5)
lines(c(1.1,1.1), CI.bounds.pca)
points(c(1.1,1.1), CI.bounds.pca, pch="_", cex=1)
mtext("Complier", 1, line=1, at=1)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

dev.off()

################################################################################

pdf("fig3-right.pdf", width=2.5, height=4.5)
par(mfrow=c(1,1), oma=c(1,1,1,1), mar=c(3.1,4,2.1,1.1))

point.late <- res.pooled.flip$t0[8]

sd.point.late <- sd(na.omit(res.pooled.flip$t[,8]))

CI.point.late <- c(min(1, max(0, point.late - qnorm(.975)*sd.point.late)), 
    min(1, max(0, point.late + qnorm(.975)*sd.point.late)))

plot(0, point.late, ylim=c(0,1), xlim=c(-.5,.5), pch=16,
    ylab = "Local Average Treatment Effect", xlab = "",
    main = "Causal Effect", xaxt="n")

lines(c(0,0), CI.point.late)
points(c(0,0), CI.point.late, pch="_", cex=1)
mtext("Complier", 1, line=1, at=0)
abline(h=0, lty="dotted")
abline(h=0.2, lty="dotted")
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

dev.off()


